function Resid = dyneqklein(Y)

% Defines equations of stochastic model;
% Input: X: All variables, current and lagged
% Outputs: Equation residuals

global Params;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LOAD parameters
 
beta = Params.beta; 
gam = Params.gam; 
chi = Params.chi; 
theta= Params.theta; 
mu = Params.mu; 
nu = Params.nu;
eta = Params.eta;
alpha = Params.alpha;

model = Params.model; 
menucost = Params.menucost; 
lbar = Params.lbar; 
noise = Params.noise; 

PMAT = Params.PMAT; 
T2 = Params.T2; 
Pgrid = Params.Pgrid; 
pstep = Params.pstep; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% KLEIN setup:  X = [Pdist, DP, mlag, V, C, PI, m]
% z is the exogenous stochastic processes driving money and productivity, with z_t+1 = phiz z_t + iidshock
% PhiHatNow is dist at beginning of t before firms adjust AND before money shock z_t is realized.
% PhiHatNow_t = Phi_t-1 = end of period t-1 production distribution
% State variables are: today's shocks z, PhiHatNow, price dispersion, and lagged real money  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% LOAD VARIABLES 
  nPhi = Params.nPhi;
  nV = Params.nV;
  % nsj = Params.nsj;
  nss = Params.nss;

% EXOGENOUS VARIABLES  
  znow = Y(Params.iz);       % exogenous shock process (possibly correlated)
  % znext = Y(Params.iznext);  % ShockNow and ShockNext are scalars: no further extraction needed
  zRnow = znow(1);           % Money growth shock today
  % zRnext = znext(1);         % Money growth shock tomorrow
  zAnow = znow(2);           % Common TFP shock

% ENDOGENOUS VARIABLES:
  xnow = Y(Params.ix);         % endogenous variables now
  xnext = Y(Params.ixnext);    % endogenous variables next
    
  % ENDOGENOUS STATE VARIABLES
  PhiHatNow = xnow(1:nPhi);
  PhiHatNext = xnext(1:nPhi);
  % DPnow = xnow(nPhi+1);
  DPnext = xnext(nPhi+1);
  mlagnow = xnow(nPhi+2);
  mlagnext = xnext(nPhi+2);

  % ENDOGENOUS JUMP VARIABLES:
  Vnow = xnow(nPhi+nss+1:nPhi+nss+nV);
  Vnext = xnext(nPhi+nss+1:nPhi+nss+nV);
  Cnow = xnow(nPhi+nss+nV+1);
  Cnext = xnext(nPhi+nss+nV+1);
  PInow = xnow(nPhi+nss+nV+2);
  PInext = xnext(nPhi+nss+nV+2);
  mnow = xnow(nPhi+nss+nV+3);
  %mnext = xnext(nPhi+nss+nV+3);
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% RESHAPE VECTORS TO MATRICES size (nump,1)
  [nump]=length(PMAT);
  siz=[nump,1];
  PhiHatNow = reshape(PhiHatNow,siz);
  PhiHatNext = reshape(PhiHatNext,siz);
  Vnow = reshape(Vnow,siz);
  Vnext = reshape(Vnext,siz);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CALCULATE SOME VARIABLES USING EQUATIONS OUTSIDE OF KLEIN REPRESENTATION

% Nominal interest rate: from money demand equation
  inow = 1/(1 - nu*Cnow^gam/mnow);

% Real wage: from Labor supply equation (detrended by P)
  wnow  = chi*Cnow^gam; 
  wnext = chi*Cnext^gam;
    
% Payoff today
  % Unow = sMAT.*Cnow.*(PMAT).^(-theta).*(PMAT - wnow*exp(-zAnow));   
  Unow = Cnow.*(PMAT.^(1-theta)-wnow*exp(-zAnow)*PMAT.^(-theta/(1-eta)));   

  
% Calculate the transition matrices R (THE OFFSET IS DETERMINED BY PI NOW)
  Rnow  = Rmatrix(nump, PInow, pstep);
  Rnext = Rmatrix(nump, PInext, pstep);

  if model>5
  % Calculate stochastic strategies logitprobMAT now and next
    logitpMATnow = logitprob(Vnow, noise, wnow);
    logitpMATnext = logitprob(Vnext, noise, wnext);
  % Calculate expected payoff EV now and next
    EVnow = sum(logitpMATnow.*Vnow);
    EVnext = sum(logitpMATnext.*Vnext);

    relentropyNow = ones(nump,1)*(log(nump) + sum(logitpMATnow.*log(max(logitpMATnow,eps^.5))));
    relentropyNext = ones(nump,1)*(log(nump) + sum(logitpMATnext.*log(max(logitpMATnext,eps^.5))));

  else
  % Get optimal value and optimal price now and next
    % [EVnow, pStarNow]   = get_pstar(Vnow, Pgrid, PInow);
    % [EVnext, ~] = get_pstar(Vnext, Pgrid, PInext);
    % PmatrixNOW = Pmatrix(pStarNow, Pgrid, pstep);
    logitprobMATnow = logitprob(Vnow, noise, wnow);
    logitprobMATnext = logitprob(Vnext, noise, wnext);
    EVnow = sum(logitprobMATnow.*Vnow);
    EVnext = sum(logitprobMATnext.*Vnext);
    
    PmatrixNOW = logitprobMATnow;
    
    relentropyNow = [];
    relentropyNext = [];
  end

  % Calculate adjustment values Dnow and Dnext
  adjcostNow = cost(model, menucost, noise, relentropyNow);
  adjcostNext = cost(model, menucost, noise, relentropyNext);

  Dnow  = ones(nump,1)*EVnow - Rnow'*T2'*Vnow - wnow.*adjcostNow;   
  Dnext = ones(nump,1)*EVnext - Rnext'*T2'*Vnext -wnext.*adjcostNext; 

  
% Calculate probabilities of adjustment Lambdanow and Lambdanext, and ExpectGAINSnext    
  LambdaNow = adjustment(model, Dnow, alpha, menucost, lbar);
  LambdaNext = adjustment(model, Dnext, alpha, menucost, lbar);
 
  timecostNext = 0;
  
  CONTINVALnext = Vnext + Dnext.*LambdaNext - wnext*timecostNext;

% Calculate PhiTildeNow from PhiHatNow = Phi_t-1
  PhiTildeNow = Rnow*T2*PhiHatNow;
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CALCULATE RESIDUALS FOR JACOBIAN

% PhiHatNext = dist of prices at time of production NOW (after shocks and adjustment)
  PhiResid =-PhiHatNext + (1-LambdaNow).*PhiTildeNow + PmatrixNOW.*(ones(nump,1)*sum(LambdaNow.*PhiTildeNow)); 
  
% Price dispersion dynamics  
  deltaResid = DPnext - sum(sum(PhiHatNext.*PMAT.^(-theta))); 
  
% Lagged real money level
  mlagResid = mlagnext - mnow;   

% Value function residual
 VResid = -Vnow + Unow + beta*(Cnext/Cnow)^(-gam)*Rnext'*T2'*CONTINVALnext;  
  %VResid = -Vnow + Rnext'*Unow + beta*(Cnext/Cnow)^(-gam)*Rnext'*T2'*CONTINVALnext;
  %VResid = -Vnow + Unow + beta*(Cnext/Cnow)^(-gam)*(Rnext'*T2'*Vnext+Dnext.*LambdaNext);
  %VResid = -Vnow + Unow + beta*(Rnext'*T2'*Vnext+Dnext.*LambdaNext);

% Consumption Euler residual
  eulerResid = mnow - Cnow; %1 - inow*beta*(Cnext/Cnow)^(-gam)/PInext; %

% Aggregate price level residual
  priceResid = 1 - (sum(sum(PhiHatNext.*(PMAT.^(1-theta)))))^(1/(1-theta));
  
% Money growth process
  moneyGrowthresid = PInow*mnow/mlagnow/mu - exp(zRnow);   
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ASSEMBLE RESID VECTOR:
  Resid = [PhiResid(:); deltaResid; mlagResid; VResid(:); eulerResid; priceResid; moneyGrowthresid];
% order of residuals is unimportant, but for simplicity follow same ordering
%   for residuals as for variables and equations
   
